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Abstract 

A  numerical  investigation  of  the  coupled  electrical  conduction  and  mass  diffusion  in  the  cathodic  GDL  of  a  PEMFC  is  performed  using  2D 
simulations.  The  current  density  on  the  GDL/catalyst  layer  interface,  which  constitutes  one  of  the  boundary  conditions  for  the  GDL  domain 
and  reflects  the  activation  overpotential  in  the  catalyst  layer  and  the  ohmic  loss  in  the  membrane,  is  solved  iteratively  using  a  novel  numerical 
algorithm.  A  parametric  study  is  performed  to  investigate  the  effects  on  current  density  distribution  of  various  operating  conditions  such  as  oxygen 
concentration  and  membrane  resistance,  and  of  design  factors  such  as  GDL  geometry,  anisotropic  transport  properties,  and  deformation  under 
the  land  area  due  to  compression.  The  results  show  that  the  current  density  distribution  under  the  land  area  can  be  dominated  by  either  electron 
transport  or  mass  transport,  depending  on  the  operating  regime.  The  analysis  of  the  in-plane  current  density  gradients  shows  the  contributions  due 
to  electrical  conduction,  oxygen  diffusion  and  membrane  resistance  in  an  explicit  form.  The  analysis  also  provides  guidance  for  the  scaling  of  the 
coupled  transport  problem. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

In  most  planar  proton  exchange  membrane  fuel  cell  (PEMFC) 
designs  the  membrane  electrode  assembly  (MEA)  is  sandwiched 
between  bipolar  plates.  The  interface  between  the  bipolar  plates 
and  the  MEA,  or  ‘land’,  is  where  electric  current  flows  in/out  of 
the  MEA.  The  impact  of  the  land  area  on  fuel  cell  operation  has 
been  investigated  numerically  by  Kulikovsky  et  al.  [1],  Meng 
and  Wang  [2],  Nguyen  et  al.  [3],  Ziegler  et  al.  [4]  and  recently 
Sun  et  al.  [5]  amongst  others.  In  these  2D  or  3D  computational 
simulations,  higher  current  density  distributions  under  the  land 
area  were  reported  for  typical  operating  conditions  and  for  a 
conventional  channel/land  geometry.  The  shift  in  current  distri¬ 
bution  patterns,  and  in  particular  of  the  maximum  current  density 
location  from  the  land  to  the  channel  area  under  some  operating 
conditions  highlights  the  trade-offs  between  electrical  and  dif¬ 
fusive  transport.  The  non-uniform  current  density  distribution 
under  the  land  associated  with  a  specific  land/channel  geometry 
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is  expected  to  have  a  significant  impact  on  the  transport  phenom¬ 
ena  in  the  lateral  direction  within  the  MEA.  Indeed  the  gradients 
in  the  lateral  direction  may  be  the  largest  in  the  entire  fuel  cell. 
From  a  design  standpoint,  understanding  the  effects  of  the  land 
is  therefore  crucial. 

Many  of  the  contributions  made  in  the  last  few  years  to  com¬ 
putational  analysis  of  PEM  fuel  cells  have  been  discussed  in  a 
thorough  review  [6]  and  notable  earlier  work  includes  Berning 
et  al.  [7],  Mazumder  and  Cole  [8],  Lee  et  al.  [9],  Meng  and 
Wang  [2],  and  Nguyen  et  al.  [3].  Some  more  recent  represen¬ 
tative  contributions  include  parallelized  three-dimensional  CFD 
modelling  [10,11,12]  parametric  computational  studies  and  val¬ 
idation  [13];  developments  addressing  the  important  aspects  of 
two-phase  transport  [14]  and  catalyst  layer  modelling  [15].  Most 
of  the  aforementioned  works  involved  numerical  solutions  of  a 
full  set  of  complex  and  coupled  transport  equations,  and  typi¬ 
cally  include  the  complete  fuel  cell  geometry.  While  allowing 
high  fidelity,  such  computations  are  intensive  and  often  imprac¬ 
tical  for  design  or  optimization  purposes.  The  aim  of  this  work 
is  to  develop  a  simplified  model  that  accounts  for  the  salient  and 
dominant  processes  associated  with  coupled  electron  and  mass 
transport  inside  the  GDL.  A  2D  configuration  is  thus  considered 
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Nomenclature 

D  mass  diffusivity  (m2  s-1) 

/  F/RUT 

F  Faraday  constant,  96487  C 

i  current  density  (A  m-2) 

z'o  cathode  exchange  current  density  (A  m-2) 

Mo2  oxygen  molecular  weight  (kg  mol  - 1 ) 

P  pressure  (Pa) 

^mem  membrane  resistance  (Q) 

Ru  universal  gas  constant,  8.3 14  J  mol-1  K-1 
T  temperature  (K) 

x  X-coordinate  (in-plane) 

xo2  oxygen  mole  fraction 

y  y-coordinate  (through-plane) 

Greek  letters 

ac  cathode  transfer  coefficient 

rj  activation  overpotential  (V) 

p  mass  density  of  mixture  (kg  m-3) 

a  electrical  conductivity  (Q ~ 1  m- 1 ) 

0  electrical  potential  (V) 

Subscripts 

BC  on  the  GDL/catalyst  layer  interface 

ref  reference  value 

vv  in-plane  direction  of  the  GDL 

yy  through-plane  direction  of  the  GDL 

Superscripts 

N  current  iteration  step 

*  on  the  anode  bipolar  plate 


for  the  GDL  taking  into  account  the  influences  of  the  catalyst 
layer  and  the  membrane  on  one  of  the  GDL  boundaries.  The 
main  objective  is  to  develop  a  model  that  adequately  describes 
the  coupling  and  that  is  suitable  for  design  and  optimization 
purposes,  where  rapid  turnaround  is  required. 

2.  Mathematical  formulation  and  numerical  method 
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Fig.  1.  Computational  domain  of  the  2D  GDL  problem. 


potential  and  zero  oxygen  flux  on  the  land,  (iii)  prescribed  oxy¬ 
gen  mole  fraction  and  zero  current  on  the  gas  channel  boundary, 
and  (iv)  prescribed  current  and  oxygen  flux  on  the  GDL/catalyst 
layer  interface: 


30  i 

3  y  CTyy 


dxo2  Mo2  l 
3 y  pD  4 F 


Boundary  condition  (iv)  requires  information  of  local  current 
density  distribution  that  satisfies 


0*  —  0bc  +  0  +  iR  (5) 

where  0  represents  the  cell  potential  at  the  anode  side  of  the 
bipolar  plate,  which  is  the  sum  of  the  ohmic  loss  in  the  GDL, 
activation  overpotential  loss  in  the  cathode  catalyst  layer,  and 
ohmic  loss  across  the  membrane.  All  other  losses  on  the  anode 
side  of  the  MEA  are  assumed  negligible.  Eq.  (5),  hereafter  called 
‘potential  balance  equation’,  is  a  boundary  condition  required 
to  close  the  model.  The  purpose  of  introducing  such  a  potential 
is  that  one  can  conveniently  include  the  effects  of  the  catalyst 
layer  and  the  membrane  into  the  solution  procedure  even  though 
only  transport  in  the  GDL  is  directly  solved.  A  Butler- Volmer 
type  of  expression  is  adopted  for  the  activation  overpotential: 


We  consider  the  transport  processes  in  the  half  GDL  domain, 
cf.  Fig.  1,  which  is  bounded  by  the  GDL/catalyst  layer  interface 
at  the  bottom,  and  by  a  half-land  and  a  half-channel  at  the  top. 
The  governing  equations  that  describe  the  transport  of  electron 
and  oxygen  are,  respectively 


d2<p  d2<p 

axxdx2  +ayydy2 

(1) 

d2xo2  d2xo2  _ 

"  dx2  +  yy  dy2 

(2) 

where  subscript  xv  denotes  the  value  in  the  in-plane  direction  and 
subscript  yy  denotes  the  value  in  the  through-plane.  The  bound¬ 
ary  conditions  for  the  governing  equations  are:  (i)  symmetry 
condition  for  the  side  boundaries  of  the  domain,  (ii)  prescribed 


It  is  noted  that  the  two  transport  processes  considered  in  the 
problem  are  only  coupled  at  the  GDL/catalyst  layer  interface. 

The  potential  loss  0  is  not  known  a  priori,  and  is  assumed  to 
be  uniform  on  the  anode  bipolar  plate.  An  iterative  method  for 
calculating  this  potential  loss  using  an  updated  current  density 
profile  based  on  local  oxygen  concentration  is  developed  in  this 
study.  For  each  iteration,  we  compute  a  new  0*  based  on  the 
results  from  previous  iteration: 
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Fig.  2.  Flow  chart  of  the  solution  procedure. 


With  this  0*  the  current  density  distribution  for  the  next  iteration 
is  then  calculated  by  substituting  (7)  into  (5): 

^  ~  ~  f-  (8) 

R 

The  advantage  of  evaluating  0*  using  (7)  is  that  the  total  current 
is  conserved  at  each  iteration.  This  can  be  shown  by  summing 
(8)  and  substituting  (7)  into  the  summation: 

■N+ 1  _  ^BC  hN 

^  R  ^  R  ^  R 


^  R  ^  R  ^  R 


A  flow  chart  of  the  solution  procedure  is  shown  in  Fig.  2.  A  com¬ 
mercial  CFD  software,  CFD-ACE+,  is  used  to  solve  the  problem 
in  conjunction  with  custom  written  user-defined  subroutines. 

3.  Results  and  discussion 

3.1.  Baseline  case 

Table  1  lists  the  values  of  the  parameters  and  dimensions  used 
for  the  parametric  study.  The  numerical  values  of  ac  and  io  are 
obtained  by  fitting  a  polarization  curve  with  a  comprehensive 
unit  cell  model,  cf.  Mazumder  and  Cole  [8].  The  total  width  of 
the  gas  channel  and  land  is  1  mm;  with  symmetry  conditions  in 
the  computation,  only  half  of  the  gas  channel  and  half  of  the 
land  need  to  be  considered.  The  effective  diffusivity  of  oxygen 
should  account  for  the  porosity  of  the  GDL  and  gas  concentra¬ 
tion,  but  for  simplicity  a  constant  value  is  used  in  the  calculation. 
The  GDL  electrical  conductivity  used  is  close  to  the  through- 
plane  conductivity  of  common  carbon  paper,  which  is  usually 
one  order  of  magnitude  lower  than  the  in-plane  value.  A  slight 
variation  of  the  membrane  resistance  in  the  in-plane  direction 
is  added  to  reflect  the  membrane  resistance  variation  due  to  the 


Table  1 

Parameter  values  used  for  the  baseline  case 


Variable 

Symbol 

Value 

Unit 

GDL  thickness 

L 

250 

[xm 

Half-channel  width 

0.5 

mm 

Half-land  width 

Wi 

0.5 

mm 

Effective  oxygen  diffusivity 

D02,eff 

in 

1 

O 

H 

X 

H 

2  -1 
m  s 

GDL  electrical  conductivity 

0 

200 

Membrane  resistance 

d0 

Amem 

in 

1 

O 

r" H 

X 

H 

Q  m2 

Pressure 

P 

3 

bar 

Temperature 

T 

343 

K 

Reference  oxygen  mole  fraction 

■*02,  ref 

0.209 

1 

Transfer  coefficient 

Oi  c 

1.15 

1 

Exchange  current  density 

io 

1.116 

Am-2 

Oxygen  concentration  exponent 

y 

1 

1 

conditions  in  the  anode  side  of  the  membrane: 

n  (  2  X  —  W  \ 

Rmem(x)  =  R°mem  (  1 - —0.15  J  (10) 

It  should  be  noted  that  most  of  the  values  listed  in  Table  1  were 
extracted  from  comprehensive  3D  simulations  of  a  complete  fuel 
cell  and  the  baseline  values  are  determined  from  the  validation  of 
these  simulations  with  experimental  data  [16].  In  the  parametric 
simulations  presented  here  the  average  current  density  is  set  to 
10,000  A  m“2. 

3.2.  Effects  of  land/gas  channel  width  ratio 

Fig.  3  shows  the  predicted  current  density  on  the  GDL/ 
catalyst  layer  interface  for  different  land:channel  ratios  (‘L:C’ 
in  the  figure).  The  baseline  case  has  an  L:C  ratio  of  5:5.  For  a 
narrower  land,  e.g.  L:C  =  3:7,  the  higher  current  density  occurs 
near  the  land  side  of  the  GDL.  As  the  land  width  is  increased,  the 
current  density  maximum  shifts  to  the  gas  channel  side  of  the 
GDL.  Since  electric  current  tends  to  flow  in  the  direction  of  least 
resistance,  the  high  current  density  region  on  the  GDL/catalyst 
layer  indicates  that  for  a  narrow  land  the  transport  of  electrons 
is  dominated  by  electrical  conduction  in  the  GDL,  while  for  the 
wider  land  case  transport  of  oxygen  becomes  dominant.  Fig.  4 
shows  the  potential  loss  of  the  cell  for  different  L:C  ratios.  The 
total  potential  loss,  which  consists  of  the  losses  in  the  GDL, 
catalyst  layer  and  membrane,  decreases  with  land  width  at  first, 
reaches  a  minimum  near  L:C  =  5:5  and  then  increases  slightly  for 
wider  land.  The  break-down  of  the  potential  losses  (note  differ- 
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Fig.  3.  Current  density  profiles  for  different  land:channel  ratios. 


X 

Fig.  5.  Current  density  profiles  for  different  oxygen  mole  fractions  in  the  gas 
channel. 


ent  scales  in  the  figure)  shows  that  the  decrease  of  total  potential 
loss  from  L:C  =  3:7  to  5:5  is  mainly  due  to  the  decrease  of  the 
loss  in  the  GDL.  Further  increase  of  the  land  width  does  not 
reduce  this  loss  significantly,  rather  the  loss  due  to  the  overpo¬ 
tential  in  the  catalyst  layer  is  increased  due  to  mass  diffusion 
constraints  caused  by  the  narrow  channel  opening. 

3.3.  Effects  of  oxygen  concentration  in  the  gas  channel 

Fig.  5  shows  the  current  density  predictions  for  different  oxy¬ 
gen  mole  fraction  in  the  gas  channel.  As  oxygen  concentration 
is  decreased,  the  diffusion  of  oxygen  becomes  more  dominant 
in  the  coupling  of  the  two  transport  processes  in  the  GDL,  and 
hence  the  high  current  density  region  moves  towards  the  gas 
channel  side  of  the  GDL  and  the  difference  of  current  density 
becomes  larger.  Such  shifts  have  been  observed  earlier  by  Meng 
and  Wang  [2]  and  Nguyen  et  al.  [3],  and  more  recently  by  Sun 
et  al.  [15].  For  a  typical  fuel  cell,  low  oxygen  concentration 


occurs  downstream  of  the  gas  channel.  Therefore  one  can  expect 
a  similar  shift  of  current  density  and  increase  of  current  density 
gradient  from  the  channel  inlet  to  the  exit  in  a  fuel  cell.  Fig.  6 
shows  the  total  potential  loss  for  the  cases  shown  in  Fig.  5.  The 
potential  loss  increases  as  oxygen  concentration  is  decreased  in 
the  gas  channel. 

3.4.  Effects  of  electrical  conductivity 

Fig.  7  shows  the  current  density  predictions  for  different  GDL 
electrical  conductivities.  Isotropic  conduction  is  assumed  for 
these  cases.  The  effect  on  the  current  density  distribution  is  sig¬ 
nificant  for  conductivity  values  lower  than  the  baseline  value. 
For  conductivity  higher  than  the  baseline  value  the  current  den¬ 
sity  is  not  sensitive  to  the  conductivity  because  the  loss  in  the 
GDL  is  negligible,  cf.  Fig.  8,  thus  the  current  density  distribution 
is  mainly  determined  by  the  transport  of  oxygen. 


Fig.  6.  Potential  losses  for  different  oxygen  mole  fractions  in  the  gas  channel. 
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Fig.  7.  Current  density  profiles  for  different  GDL  electrical  conductivities. 

3.5.  Effects  of  membrane  resistance 

The  conduction  of  protons  in  the  membrane  is  not  directly 
modeled  here,  but  appears  instead  as  a  potential  loss  term  in  the 
potential  balance  Eq.  (5).  This  implies  that  the  current  flow  in  the 
membrane  is  perpendicular  to  the  GDL/catalyst  layer  interface. 
In  this  sense  the  effect  of  membrane  resistance  can  be  consid¬ 
ered  as  an  additional  through-plane  resistance  to  the  electrical 
conduction  in  the  GDL.  Fig.  9  shows  the  predictions  of  current 
density  distribution  for  different  membrane  resistance  values, 
marked  as  the  ratio  with  respect  to  the  baseline  value  Ro  in  the 
figure.  One  can  see  that  increasing  the  membrane  resistance  has 
a  similar  effect  as  decreasing  the  GDL  electrical  conductivity, 
cf.  Fig.  7.  A  similar  trend  should  be  expected,  though  this  was 
not  ascertained  in  the  present  study,  for  the  contact  resistance  at 
the  land/GDL  interface. 

Fig.  10  shows  the  potential  loss  increases  linearly  with  mem¬ 
brane  resistance  values.  The  implication  of  the  parametric  study 


1/sigmaY 

Fig.  8.  Potential  losses  for  different  GDL  electrical  conductivities. 


Fig.  9.  Current  density  profiles  for  different  membrane  conductivities. 


with  membrane  resistance  is  that  for  a  fuel  cell  that  operates 
under  low  humidity  conditions,  such  as  air-breathing  fuel  cell 
or  a  cell  with  low  anode  humidification,  the  current  density 
distribution  in  the  GDL  may  be  quite  different  from  that  in  a 
well-humidified  fuel  cell. 

3.6.  Effects  of  compression 

In  assembling  the  bipolar  plate  and  the  MEA  in  a  fuel  cell, 
a  compression  force  is  always  required  to  ensure  proper  seal¬ 
ing  and  to  reduce  contact  resistance  between  the  bipolar  plate 
and  the  GDL.  The  GDL  is  expected  to  undergo  some  deforma¬ 
tion  under  such  compression  as  it  is  the  least  stiff  component 
in  the  MEA.  The  impact  of  the  compression  on  transport  pro¬ 
cesses  in  the  GDL,  in  particular  the  region  underneath  the  land, 
is  likely  two-fold:  the  through-plane  conductivity  of  the  GDL 
would  increase  as  new  conducting  pathways  form  between  the 
fibres,  and  the  through-plane  diffusivity  of  oxygen  would  be 
reduced  as  the  porosity  is  reduced.  Based  on  the  observations  in 


Fig.  10.  Potential  losses  for  different  membrane  conductivities. 
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Fig.  11.  Current  density  profiles  for  different  GDL  conductivities  and  diffusiv- 
ities. 


the  aforementioned  parametric  studies,  both  changes  will  pro¬ 
mote  electrical  conduction  in  the  GDL  and  oxygen  diffusion 
will  become  dominant.  Fig.  1 1  shows  the  numerical  results  of 
current  density  distribution  with  properties  changed  only  for 
the  portion  of  GDL  under  the  land  to  simulate  the  effects  of 
compression.  The  bold  line  (a)  is  the  baseline  case,  which  uses 
isotropic  properties.  When  the  diffusivity  is  reduced  to  half  and 
the  conductivity  doubled  (b)  for  an  isotropic  case  and  (c)  for  an 
anisotropic  case,  departure  from  the  baseline  is  observed.  Fur¬ 
ther  reducing  the  diffusivity  and  increasing  the  conductivity  by 
an  order  of  magnitude  (d  and  e)  shows  similar  trends  but  more 
pronounced  effects.  Fig.  12  compares  the  oxygen  mole  fraction 
distribution  and  the  current  flow  in  the  GDL  for  the  baseline  case 
and  case  (e)  in  Fig.  1 1 ,  which  has  a  low  through-plane  diffusivity 
(1/10  of  the  baseline  value)  and  a  high  through-plane  conductiv¬ 
ity  (10  times  of  the  baseline  value)  in  the  under- the-land  portion 
of  the  GDL.  For  case  (e)  the  oxygen  diffusion  into  the  com- 
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Fig.  12.  Oxygen  mole  fraction  in  the  GDL  for  different  compression  conditions 
(top:  baseline,  bottom:  case  (e)  in  Fig.  1 1). 


pressed  region  is  hindered  by  the  low  diffusivity,  hence  the  low 
current  density  in  this  region.  Due  to  the  relatively  high  conduc¬ 
tivity  in  the  region  under  the  land,  most  of  the  current  under  the 
channel  flows  right  across  the  border  between  the  compressed 
region  and  the  uncompressed  region. 

3.7.  Analysis 


Comprehensive  3D  CFD-based  simulations  of  the  PEMFC 
have  the  advantage  of  providing  detailed  information  of  the 
complex  transport  phenomena  involved  in  the  system.  How¬ 
ever,  using  comprehensive  CFD  simulation  tools  to  gain  physical 
insights  and  scaling  of  the  problem  is  computationally  expen¬ 
sive.  The  simplified  model  of  the  present  study  provides  a  rapid 
method  for  analyzing  the  coupled  transport  in  the  ME  A.  In  the 
present  model,  the  transport  of  gas  species  and  electron  in  the 
GDL  are  both  described  by  two  Laplace  equations  with  similar 
boundary  conditions.  The  coupling  of  both  transport  processes 
occurs  only  on  one  of  the  domain  boundaries,  i.e.  on  the  inter¬ 
face  of  the  cathodic  catalyst  layer  and  the  GDL.  The  effects  of 
the  catalyst  layer  and  the  membrane  are  captured  through  the 
potential  balance  Eq.  (5),  which  essentially  is  a  boundary  condi¬ 
tion  for  electrical  potential  when  the  entire  unit  cell  (MEA  and 
bipolar  plates)  is  considered.  In  the  model,  the  potential  balance 
equation  is  used  to  evaluate  the  current  density  profile  on  the 
boundary  where  coupling  of  gas  species  and  electron  occurs. 
In  fuel  cell  operation,  one  is  concerned  with  the  distribution  of 
current  in  the  MEA.  A  current  density  distribution  of  minimum 
spatial  variation  in  general  yields  an  optimal  cell  performance 
and  implies  maximum  durability.  The  gradient  of  the  through- 
plane  current  density  in  the  in-plane  direction  on  the  catalyst 
layer-GDL  interface  can  be  obtained  from  the  potential  balance 
Eq.  (5)  by  taking  the  derivative  of  the  potential: 


90*  _  90b c  dr]  d(iR ) 

dx  dx  dx  dx 


(ID 


By  applying  the  chain  rule  to  the  activation  overpotential  and  the 
ohmic  terms,  and  with  the  assumption  of  a  uniform  cell  potential 
0  ,  we  have 
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(12) 


With  (6)  for  the  activation  overpotential  and  rearranging  (12), 
we  have  the  gradient  of  current  density  as 


di  —  1  /  90bc  .  dR  —  1  9vo2 

dx  R  +  — \  dx  dx  xo2acf  dx 

The  three  terms  in  the  bracket  of  (13)  represent  respectively 
the  contributions  due  to  gradient  of  electrical  conduction  in  the 
GDL,  variation  of  protonic  resistance  in  the  membrane,  and  gra¬ 
dient  of  oxygen  distribution  in  the  GDL.  Fig.  1 3  shows  the  profile 
of  these  terms  and  the  profile  of  predicted  current  density  for  the 
baseline  case.  One  can  see  that  the  current  density  reaches  a 
maximum  or  minimum  at  the  location  where  the  summation  of 
the  three  terms  in  (13)  equals  zero.  The  implication  of  the  analy¬ 
sis  on  the  current  density  gradient  is  that  in  order  to  minimize  the 
gradient  of  current  density  in  the  MEA,  one  must  (a)  minimize 
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Fig.  13.  Profiles  of  current  density  and  gradient  of  current  density  on  the  inter¬ 
face  of  the  GDL  and  catalyst  layer. 

the  sum  of  contributions  due  to  the  three  major  sources  of  non¬ 
uniformity  in  (13),  and  (b)  increase  the  membrane  resistance, 
which  is  undesirable.  We  note  that  the  first  and  third  terms  in 
the  bracket  of  (13)  are  of  opposite  sign  under  normal  conditions. 
When  the  land  area  is  wide,  or  for  a  GDL  material  with  low  in¬ 
plane  conductivity,  the  contribution  due  to  electrical  conduction 
in  the  GDL  is  relatively  small  compared  to  that  due  to  mass  dif¬ 
fusion  in  the  GDL.  In  this  case  the  current  gradient  is  dominated 
by  the  oxygen  distribution  on  the  catalyst  layer  surface.  On  the 
other  hand,  for  a  large  channel  opening,  the  effects  due  to  oxygen 
distribution  become  less  significant  and  the  electrical  conduc¬ 
tion  may  become  dominant.  The  balance  between  the  electrical 
conduction  and  mass  diffusion  in  the  GDL  is  also  affected  by  the 
distribution  of  membrane  resistance,  which  is  closely  related  to 
current  density  current  density  distribution  and  local  humidity. 

4.  Conclusions 

In  this  work  we  present  a  simplified  model  for  simulating 
coupled  electron  and  mass  transport  in  the  GDL  of  a  PEMFC. 
The  parametric  study  performed  using  this  model  shows  that 
current  density  distribution  at  the  GDL/catalyst  interface  can 
be  dictated  by  either  electron  transport  or  mass  transport  in 
the  GDL  depending  on  parameters  including  channel  geome¬ 
try  (land  versus  channel  width),  oxygen  concentration  in  the 
gas  channel,  transport  properties  (conductivity  and  diffusiv- 
ity),  membrane  resistance,  etc.  The  model  is  also  capable  in 
predicting  the  effects  due  to  compression,  provided  a  realistic 
correlation  between  compression  force  and  anisotropic  transport 
properties  is  available.  An  analysis  based  on  examining  the  gra¬ 
dients  on  the  catalyst  layer  surface  has  shown  that  the  current 
gradient  is  dependent  on  the  gradients  of  electrical  potential  in 
the  GDL,  oxygen  concentration  and  water  content  in  the  mem¬ 
brane.  The  approach  presented  here  provides  an  effective  method 


for  assessing  design  trade-offs  and  for  optimization  of  channel 
dimensions.  Ongoing  and  future  work  [17]  will  focus  on  using 
this  model  to  provide  design  guidelines. 
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